function [f,x,xp,y,yp,symparams,Phi,h1,nomSDF_cup,heteroShock,r_cu] = DSGEmodel(unitFree)
% This function reports the equations for the New Keynesian model with Calvo prices
% in f, and it reports the state vector (x and xp) and the control vector (y and yp)

% IMPORTANT: lagged controls leaded one period forward must be written as
% the current controls. That is c_cu and not c_ba1p

%Define parameters
syms BETTA B CHI CHI0 THETA DELTA ALFA PHI PHIzero ZI KAPAw ETA RHOR PHIpai PHIpai_1 ...
     PHIy PHIy_1 PHIc PHIc_1 PHIl NU U0 U0d...
     OMEGAz OMEGAd OMEGAn OMEGAp OMEGAa ...
     RHOz  RHOd   RHOn  RHOp RHOa ...
     Kss OUTPUTss PAIss MUZss Css AA lss Wss;
Rss = sym('Rss');
symparams=[BETTA,B,CHI,CHI0,THETA,DELTA,ALFA,PHI,PHIzero,ZI,KAPAw,ETA,...
     RHOR,PHIpai,PHIpai_1,PHIy,PHIy_1,PHIc,PHIc_1,PHIl,NU,U0,U0d,...
     OMEGAz,OMEGAd,OMEGAn,OMEGAp,OMEGAa,...
     RHOz ,RHOd   RHOn  ,RHOp,RHOa, ...
     Kss OUTPUTss PAIss MUZss Css AA lss,Rss,Wss];

%Define variables 
syms c_cu   pai_cu   evf_cu  ...
     c_cup  pai_cup  evf_cup ...
     c_ba1  muz_cu   d_cu  n_cu  paistar_cu  a_cu...
     c_ba1p muz_cup  d_cup n_cup paistar_cup a_cup;
 

%EQ7: Agg. ressource constraint
%f0 = - output_cu*(1-ZI/2*(pai_cu /PAIss^NU-1)^2) + (c_cu  + DELTA*Kss);
output_cu  = (c_cu  + DELTA*Kss)/(1-ZI/2*(pai_cu /PAIss^NU-1)^2);
output_cup = (c_cup + DELTA*Kss)/(1-ZI/2*(pai_cup/PAIss^NU-1)^2);

% EQ 8: The aggregated production function to get output. 
%output_cu  = a_cu *Kss^THETA*l_cu ^(1-THETA);  
l_cu  = (output_cu /(a_cu *Kss^THETA))^(1/(1-THETA));
l_cup = (output_cup/(a_cup*Kss^THETA))^(1/(1-THETA));

% EQ 2: Household's First-order condition for labour
%f1= -n_cu*PHIzero*(1-l_cu)^(-1/PHI) + ((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(-CHI)*w_cu/Css^CHI0;
w_cu = KAPAw*Wss+(1-KAPAw)*(n_cu*PHIzero*(1-l_cu)^(-1/PHI)/((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(-CHI)*Css^CHI0);

% EQ 4: Firms FOC for labour
%mc_cu  = w_cu/((1-THETA)*a_cu *Kss^THETA*l_cu ^-THETA);
mc_cu  = w_cu/((1-THETA)*a_cu *Kss^THETA*l_cu ^-THETA);

% Taylor rule
r_cu = Rss*exp((PHIpai*(log(pai_cu)-(log(PAIss) + log(paistar_cu)) )...
         + PHIpai_1*(log(pai_cup)-(log(PAIss) + log(paistar_cup))) ...
         + PHIy*log(output_cu/OUTPUTss)...
         + PHIy_1*log(output_cup/OUTPUTss)...
         + PHIc*(c_cu-c_ba1+log(muz_cu)-log(MUZss)) ...
         + PHIc_1*(c_cup-c_ba1p+log(muz_cup)-log(MUZss)) ...
         + PHIl*(log(l_cu/lss))));

% The expression for the value function (minus)
% vf_cu = - (U0 + d_cu*(1/(1-CHI)*((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(1-CHI)+U0d + n_cu*PHIzero/(1-1/PHI)*(1-l_cu )^(1-1/PHI)) ) ...
%         + BETTA*AA*evf_cu ^(1/(1-ALFA));
vf_cup = - (U0 + d_cup*(1/(1-CHI)*((c_cup-B*c_cu/muz_cup)/Css^CHI0)^(1-CHI)+U0d + n_cup*PHIzero/(1-1/PHI)*(1-l_cup)^(1-1/PHI)) ) + BETTA*AA*evf_cup^(1/(1-ALFA));    

% The real stochastic discount factor
mReal_cup = BETTA*d_cup/d_cu*...
    (AA*evf_cu^(1/(1-ALFA))/(vf_cup*muz_cup^((1-CHI)*(1-CHI0)) ))^ALFA...   
    *(c_cup-B*c_cu/muz_cup)^(-CHI)/((c_cu-B*c_ba1/muz_cu)^(-CHI))*muz_cup^(-CHI*(1-CHI0)-CHI0);  

%% The equations in f 
% EQ 3: Euler-equation for the one-period interest rate
if unitFree == 1
    f1 = -1 + mReal_cup/pai_cup*r_cu;
else
    f1 = -1 + mReal_cup/pai_cup*r_cu;
end

% EQ 5: First-order condition for prices
if unitFree == 1
    f2 = ((1-ETA)+ZI*mReal_cup*(pai_cup/PAIss^NU-1)*pai_cup/PAIss^NU*output_cup/output_cu*muz_cup ...
        - ZI*(pai_cu/PAIss^NU-1)*pai_cu/PAIss^NU)/(ETA*mc_cu) + 1;
else
    f2 = (1-ETA)+ZI*mReal_cup*(pai_cup/PAIss^NU-1)*pai_cup/PAIss^NU*output_cup/output_cu*muz_cup ...
        + ETA*mc_cu - ZI*(pai_cu/PAIss^NU-1)*pai_cu/PAIss^NU;
end

% The expected value of the value function 
if unitFree == 1
    f3 = -1 + (vf_cup*muz_cup^((1-CHI)*(1-CHI0))/AA)^(1-ALFA)/evf_cu;
else
    f3 = -evf_cu + (vf_cup*muz_cup^((1-CHI)*(1-CHI0))/AA)^(1-ALFA);
end

%% The links
% For lagged consumption
if unitFree == 1
    f4 = -1 + c_cu/c_ba1p;
else
    f4 = -c_ba1p + c_cu;
end

%% The linear shocks
%For non-stationary technology
f5 = -log(muz_cup/MUZss) + RHOz*log(muz_cu/MUZss);

%The preference shock
f6 = -log(d_cup) + RHOd*log(d_cu);

%For labor shocks
f7 = -log(n_cup) + RHOn*log(n_cu);

%Inflation target
f8 = -log(paistar_cup) + RHOp*log(paistar_cu);

%stationary technology shocks
f9 = -log(a_cup) + RHOa*log(a_cu);

%Creating function f
f =[f1; f2; f3; f4; f5; f6; f7; f8; f9];

% Define the variables which needs to be log-transformed
y  = [c_cu    pai_cu   evf_cu   ];
yp = [c_cup   pai_cup  evf_cup  ];
x  = [c_ba1  muz_cu  d_cu  n_cu   paistar_cu   a_cu];
xp = [c_ba1p muz_cup d_cup n_cup  paistar_cup  a_cup];

% Make f a function of the logarithm of the state and control vector
f = subs(f, [x,y,xp,yp], exp([x,y,xp,yp]));

h1 = [];    %endogenous states (not lagged controls)
Phi = [];
nomSDF_cup  = subs(mReal_cup/pai_cup,[x,y,xp,yp], exp([x,y,xp,yp]));
heteroShock = [OMEGAz;OMEGAd;OMEGAn;OMEGAp;OMEGAa];

% Re-defining the controls and the states - if needed

% Extra
r_cu = subs(r_cu, [x,y,xp,yp], exp([x,y,xp,yp]));



end